setwd("C:/Users/shengwuxinxi/Desktop/file")
file1=read.table("2_5hmc_yjp_bam/DMR/DC6v6q_cc8v8q_anno.txt",head=T,sep="\t",fill=T,quote="")
file2=read.table("13968DMRs_9v9p_8v8p.txt",header = T,sep="\t")
file1$chr_ss=paste(file1$chr,file1$start,file1$stop,sep="_")

file11=file1[,1:12]
file22=file2[,1:2]
file3=merge(file11,file22,by="chr_ss")
id=as.character(file3$chr_ss)
row.names(file1)=file1$chr_ss
row.names(file2)=file2$chr_ss
file11=file1[id,]
file22=file2[id,]

sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
sel2=c("M30_M29","M26_M25","M35_M36","M28_M27")
sel1=unlist(strsplit(sel1,"_"))[seq(2,2*length(sel1),2)]
sel2=unlist(strsplit(sel2,"_"))
sel=c(paste(sel1,".counts",sep=""),paste(sel2,".counts",sep=""))

file12=data.frame(file11[,1:9],file11[,sel])
a=file12[,sel]

sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M42_M41","M44_M43","M48_M47","M50_M49","M52_M51")
sel2=c("M30_M29","M26_M25","M35_M36","M28_M27")
sel1=unlist(strsplit(sel1,"_"))[seq(2,2*length(sel1),2)]
sel2=unlist(strsplit(sel2,"_"))
sel=c(paste(sel1,".counts",sep=""),paste(sel2,".counts",sep=""))

file23=data.frame(file22[c("chr_ss",sel)])
b=file23[,sel]


file12$mean=0
file23$mean=0
for (i in 1:39){
file12[i,]$mean=mean(unlist(a[i,]),na.rm=T)
file23[i,]$mean=mean(unlist(b[i,]),na.rm=T)
}

file12$chr_ss=paste(file12$chr,file12$start,file12$stop,sep = "_")
rt=data.frame(chr_ss=file12$chr_ss,genename=file12$Gene.Name)
file23=merge(rt,file23,by="chr_ss")

a=data.frame(chr_ss=file12$chr_ss,genename=file12$Gene.Name,mean=file12$mean)
a=a[order(a$mean),]
a$group="5hmc"
a$num=1:39
c=data.frame(chr_ss=a$chr_ss,num=a$num)

b=data.frame(chr_ss=file23$chr_ss,genename=file23$genename,mean=file23$mean)
b$group="5mc"
b=merge(b,c,by="chr_ss")
cc=rbind(a,b)
ggplot(cc,aes(x=cc$num,y=cc$mean))+scale_x_continuous(breaks = cc$num,labels = cc$genename)+labs(x="genename",y="mean")+
geom_point(aes(color=group),size=1)+geom_line(aes(color=group),size=1)+
theme_bw()+theme(axis.text.x = element_text(angle = 45,hjust=1,size=10),panel.grid.minor = element_blank())
